# -*- coding:utf-8 -*-


import numpy as np



nmax = 10

#f(U)+H(U)*V=0


def evalFunction(f, U):
	"""Evaluate a function with the given vector (deprecated)"""
	evalu = np.zeros(f.shape);
	n = f.shape[1]
	m = f.shape[0]
	for i in range(m):
		for j in range(n):
			evalu[i,j] = f[i,j](U)
	return evalu

#Old Version which takes lambdas in parameter 
#def newton_raphson(f, J, U0, N, epsilon):
#	n = J.shape[1]
#	m = J.shape[0]
#	U = U0
#	for i in range(N):
#		fU = evalFunction(f,U)
#		JU = evalFunction(J,U)
#		if(np.linalg.norm(fU) < epsilon):
#			return U
#		V = np.linalg.lstsq(JU,-fU)[0]
#		U = U + V
#	return U

def newton_raphson_backtrack(f, J, U0, N, epsilon):
	"""Newton Raphson Algorithm with the backtracking optimisation"""
	U = U0
	for i in range(0,N):
		step = 1.0
		normfU = np.linalg.norm(f(U))
		if (normfU < epsilon):
			break
		V = np.linalg.lstsq(J(U), -f(U))[0]
		while(np.linalg.norm(f(U + (step*V))) > normfU): 
			step *= 2./3
		U = U + step*V
	return U



def newton_raphson_simple(f, J, U0, N, epsilon):
	"""Simple Newton Raphson Algorithm"""
	U = U0
	for i in range(0,N):
		if ((np.linalg.norm(f(U))) < epsilon):
			break
		V = np.linalg.lstsq(J(U), -f(U))[0]
		U = U + V
	return U


def newton_raphson(f ,J, U0, N, epsilon):
	"""Interface with others projects (use the backtracking method)"""
	#return newton_raphson_simple(f, J,U0, N, epsilon)
	return newton_raphson_backtrack(f, J, U0, N, epsilon)


